First-principles study of spontaneous polarization in multiferroic BiFeC>3 
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The ground-state structural and electronic properties of ferroelectric BiFeOs are calculated using 
density functional theory within the local spin-density approximation and the LSDA+U method. 
The crystal structure is computed to be rhombohedral with space group R3c, and the electronic 
structure is found to be insulating and antiferromagnetic, both in excellent agreement with available 
experiments. A large ferroelectric polarization of 90-100 /iC/cm 2 is predicted, consistent with the 
large atomic displacements in the ferroelectric phase and with recent experimental reports, but dif- 
fering by an order of magnitude from early experiments. One possible explanation is that the latter 
may have suffered from large leakage currents. However both past and contemporary measurements 
are shown to be consistent with the modern theory of polarization, suggesting that the range of re- 
ported polarizations may instead correspond to distinct switching paths in structural space. Modern 
measurements on well-characterized bulk samples are required to confirm this interpretation. 

PACS numbers: 71.20.-b, 77.80.-c 
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I. INTRODUCTION 

There has been considerable recent interest in develop- 
ing multifunctional materials in which two or more use- 
ful properties are combined in a single compound. Per- 
haps the most widely studied class of multifunctional ma- 
terials are the diluted magnetic semiconductors, where 
interaction between magnetic and electronic degrees of 
freedom allows both charge and spin to be manipu- 
lated by applied homogeneous electric fields i These and 
other related materials having spin-dependent electronic 
properties are presently being explored for spintronic 
applications. 2 However another class of materials, the 
so-called multiferroicsj^ is also of growing importance. 
Multiferroic materials have simultaneous ferromagnetic, 
ferroelectric and/or ferroelastic ordering. Coupling be- 
tween the magnetic and ferroelectric order parameters 
can lead to magnetoelectric effects, in which the magne- 
tization can be tuned by an applied electric field and vice 
versa. Relatively few multiferroics have been identified, 5 
and in those that are known, the mechanism underlying 
their ferroelectricity is often unconventional. 6-7 » The pur- 
pose of this work is to understand the unusual ferroelec- 
tric behavior in multiferroic BiFeOs, which has recently 
emerged as an especially promising magnetoelectric mul- 
tiferroic material^ 

Recently, large ferroelectric polarizations, exceed- 
ing those of prototypical ferroelectrics BaTiC>3 and 
PbTiC>3, have been reported in high quality thin films 
of BiFe03&2iiS These sizeable polarizations are con- 
sistent with the observed large atomic distortions?ii*i£ 
but apparently inconsistent with earlier studies of bulk 
BiFeOsfi 3 " a difference whose origin is currently un- 
der debate. In addition, appreciable magnetizations 
(~ 1 /ie/formula unit), increasing with decreasing film 
thickness, have been reported^ accompanied by substan- 



tial magnetoelectric coupling. 

Bulk BiFeOs has long been known to be ferroelectrioi 3 . 
with a Curie temperature of about 1100 K. The struc- 
ture of the ferroelectric phase, resolved experimentally 
using both X-ray and neutron diffraction^ii^ can be 
understood as highly-distorted perovskite with rhombo- 
hedral symmetry and space group i?3c. The primitive 
unit cell contains two formula units (10 atoms) as shown 
in Fig. ^ The i?3c symmetry permits the development 
of a spontaneous polarization along [111], and Bi, Fe, 
and O are displaced relative to one another along this 
3-fold axis. The largest relative displacements are those 
of Bi relative to O, consistent with a stereochemically- 
active Bi lone pairi The polar displacements (relative 
to cubic perovskite) are noticeably extreme when com- 
pared with those in non-lone-pair-active perovskite fer- 
roelectrics such as BaTiC>3 or KNbC>3, but are consis- 
tent with those observed in other Bi-based perovskites. 1 - 
The observed counter-rotations of neighboring octahedra 
about [111] are consistent with a tolerance factor some- 
what less than unity. The crystal structure of the para- 
electric phase has not been conclusively determined. 

Surprisingly, given the large atomic displacements rel- 
ative to the centrosymmetric cubic perovskite structure, 
and the high ferroelectric Curie temperature, early mea- 
surements on bulk single crystals^ yielded rather small 
polarizations. Teague et alm& initially reported a polar- 
ization along [111] of just 6.1 /j,C/cm 2 ; although these 
authors state that their hysteresis loops were not satu- 
rated, their measurements were supported by more recent 
studies of both BiFeOs/BaTiOs alloy filmsi 5 - and pure 
BiFeC>3 polycrystalline films;— The small values are in 
sharp contrast with recent experiments on epitaxial thin 
film samples of BiFeOs, which were found to possess large 
polarizations. The first thin film measurements 8 yielded 
values of 50-90 /iC/cm 2 on (100) aligned substrates, in- 
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FIG. 1: Structure of R3c BiFeOs. Notice the position of the 
oxygen octahedra relative to the Bi framework; in the ideal 
cubic perovskite structure the oxygen ions would occupy the 
face-centered sites. 



TABLE I: Various measured values for the polarization in 
BiFeOa, in chronological order with the oldest at the top. 



Ref. 


P (/xC/cm 2 ) 


sample type 


13 


6.1 


bulk single crystals 


F5 


2.5 


(Bi .7Bao.3)(Feo.7Tio.3)0 3 films (300 nm) 
on Nb-doped SrTiOs 


16 


2.2 


polycrystalline films (200 nm) 




50 - 90 


thin films (400 - 100 nm) 
on SrRuOa/SrTiOs 




35.7 


polycrystalline films (350 nm) 


XL 


8.9 


bulk ceramics 


M 


158 


polycrystalline films (300 nm) 



creasing to 100 /zC/cm 2 for (111) orientations, 9 Follow- 
ing that, a whole variety of different experimental values 
have been reportedi2*i2ii^ (summarized in Table P) , in- 
cluding a very recent report of giant (> 150 /zC/crri 2 ) 
polarization which is the highest value ever measured for 
a ferroelectric. 

There are several plausible explanations for the spread 
of experimental values. First, the original reports of 
small polarization might have been limited by poor sam- 
ple quality, with the large thin film values representing 
the "true" polarization for i?3c BiFeOa . A second possi- 
bility is that the small values could be correct for the i?3c 
structure, with the large values being correct for different 
structural modifications stabilized in the thin films. And 
finally, a third possibility is that large and small values 
can be explained within the modern theory of polariza- 
tion, which recognizes that polarization is in fact a lattice 
of values, rather than a single vector ,i2i22i2I In this latter 
case, the ferroelectric switching behavior in the different 
samples would have to be substantially different. In this 
work, we use first-principles density functional calcula- 
tions to examine this issue carefully. We find that the 
most natural value of the polarization is 90-100 /iC/cm 2 



along the [111] direction, consistent with recent thin film 
measurements A* However both the unexpectedly small 
early values, and the anomalously large recent values, can 
be explained within the modern theory of polarization, 
provided that a suitable switching path can be found. 

The remainder of this paper is organized as follows. 
In Section ^ we describe the ab initio methods used 
throughout this work. In Sections 1111 Al and 1111 BI we 
report the computed ground state structural and elec- 
tronic properties of BiFeC>3 and show that they are in 
good agreement with experiment. In Section IIII CI the 
ferroelectric polarization is calculated using the modern 
theory of polarization and compared with a simple esti- 
mate. In section IIVI we discuss how the dependence of 
the polarization on the switching path is a possible ex- 
planation for the large spread of reported experimental 
values. The intriguing magnetic properties are discussed 
in a separate studyi^ 



II. METHOD 

To calculate the structure, polarization and Born ef- 
fective charges of BiFeOs, we use density functional the- 
ory (DFT) within the local spin-density approximation 
(LSDA)22ik and the LSDA+U method 2 ^ as implemented 
in the Vienna ab initio Simulation Package (VASP)^§i2£ 
All results were obtained using the projector-augmented 
plane-wave (PAW) method 28-29 by explicitly treating 
15 valence electrons for Bi (5e? 10 6s 2 6p 3 ), 14 for Fe 
(3p 6 3c? 6 4s 2 ), and 6 for oxygen (2s 2 2p 4 ). Our scalar- 
relativistic calculations do not include spin-orbit correc- 
tions. All structural relaxations are performed within 
the LSDA. The ions are steadily relaxed toward equilib- 
rium until the Hcllmann-Feynman forces are less than 
10~ 3 eV/A. Brillouin zone integrations are performed 
with a Gaussian broadening^ of 0.1 eV during all relax- 
ations. These calculations are performed with a 3x3x3 
Monkhorst-Pack k-point mesh^i centered at T and a 
500 eV plane-wave cutoff, both of which result in good 
convergence of the computed ground state properties. 

Since it is well known that the LSDA often underesti- 
mates the size of the band gap in systems with strongly 
localized d orbitals, and even predicts metallic behav- 
ior for materials that are known to be insulators^ we 
also calculate the electronic structure (using the struc- 
tural parameters obtained within the LSDA) within the 
LSDA+U method^ In the LSDA+U framework the 
strong Coulomb repulsion between localized d states is 
treated by adding a Hubbard-like term to the effective 
potential, leading to an improved description of corre- 
lation effects in transition metal oxides. The LSDA+U 
method requires two parameters, the Hubbard param- 
eter U and the exchange interaction J. Since there is 
no unique way of including a Hubbard term within the 
DFT-framework, several different approaches exist which 
all give similar results. In this work we use the approach 
described by Dudarev et al£& where only an effective 
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Hubbard parameter U c s — U — J enters the Hamilto- 
nian. The magnitude of U c g is varied between eV and 
7 eV for the Fe d states (the standard LSDA result corre- 
sponds to U e s — eV). These calculations are performed 
with a r centered 5x5x5 Monkhorst-Pack k-point mesh 
and a slightly lower plane-wave cutoff of 400 eV. The 
tetrahedron method^* is used for Brillouin zone integra- 
tions. 

Rock salt, or G-type, antiferromagnetic (AFM) order 
is assumed for all calculations, as well as a homogeneous 
and collinear spin arrangement. This assumption is well 
justified, since in practice BiFe03 is observed to be nearly 
G-type AFM. However, experiments also report a long 
wavelength spiral spin structure 3 ^ and possibly a small 
out-of-plane canting due to weak ferromagnetism, 36 Since 
the noncollinearity is quite minimal, the simplification to 
a collinear magnetic structure is acceptable. Detailed cal- 
culations of the effects of noncollinearity and spin-orbit 
coupling on the magnetic properties of BiFe03 appear 
separately^ 

The electronic contribution to the polarization is calcu- 
lated as a Berry phase using the method first developed 
by King-Smith and Vanderbihi^ (see also Ref.Hj, the 
so-called "modern" theory of polarization. In this ap- 
proach, the total polarization P for a given crystalline 
geometry can be calculated as the sum of ionic and elec- 
tronic contributions. The ionic contribution is obtained 
by summing the product of the position of each ion in the 
unit cell (with a given choice of basis vectors) with the 
nominal charge of its rigid core. The electronic contribu- 
tion to P is determined by evaluating the phase of the 
product of overlaps between cell-periodic Bloch functions 
along a densely-sampled string of neighboring points in 
k-space. Here we use 4 symmetrized strings consisting of 
15 k points to obtain the electronic contribution to the 
polarization, which is calculated separately for each spin 
channel; the total polarization is then the sum of the two 
spin contributions and the ionic contribution. 



III. RESULTS AND DISCUSSION 

A. Structure 

In Table ITT1 we report structural parameters obtained 
by relaxing the cell volume f2, rhombohedral angle a, 
and atomic positions within the R3c space group, which 
are all in good agreement with the experimental values 
of Kubel and Schmidt The lattice constant is under- 
estimated by roughly 3%, a large but not atypical con- 
sequence of the LSDA. The internal atomic coordinates 
are reproduced very well. Both the experimental and the 
theoretical value are very close to 60° which would corre- 
spond to perfectly cubic lattice vectors. The calculated 
rhombohedral angle is slightly larger than experiment; 
since it is found to increase with decreasing volume, we 
attribute this to the underestimated volume of the LSDA 
ground state. Our calculations of the total energy for 



other selected symmetries, including PAmm, i?3c, and 
R3m, corroborate R3c as the ground state. 

TABLE II: Calculated and measured structural parameters of 
BiFeOs in space group R3c (point group Cz v ); the Wyckoff 
positions 2a and 6b are referenced to the rhombohedral sys- 
tem and are Bi(x,x,x), Fe(x,x,x), and 0(x,y,z). Also in- 
cluded are the lattice constant a r h (rhombohedral unit cell), 
the rhombohedral angle a, and the unit cell volume Q. 







LSDA 


Exp. (Ref.^2) 


Bi (2a) 


X 








Fe (2a) 


X 


0.231 


0.221 


O (6b) 


x 


0.542 


0.538 




V 


0.943 


0.933 




z 


0.408 


0.395 


a rb (A) 




5.46 


5.63 


a(°) 




60.36 


59.35 


n (A 3 ) 




115.98 


124.60 



In i?3c, the Bi site is strongly distorted such that only 
6 of the 12 oxygens surrounding Bi can still be consid- 
ered nearest neighbors; three co-planar oxygens lie above 
Bi along [111] at 2.30 A, and three sit below at 2.41 A. 
Likewise, the Fe site is displaced relative to the center of 
its surrounding octahedron, which is also distorted, with 
3 oxygen neighbors at 1.92 A and 3 others at 2.07 A. 
The O-Fe-0 bond angle in this system, which would be 
an ideal 180° in a cubic perovskite structure, buckles to a 
value of 165° in i?3c. Indeed i?3c BiFeOs is structurally 
similar to ferroelectric LiNb03, which also has i?3c sym- 
metry, off-centering of the Li and Nb cations, and large 
rotations of the oxygen octahedra resulting in a 6-fold 
coordinated site for the large cationi^LiS 

The structural relationship between cubic perovskite 
and i?3c can be understood with just two rather sim- 
ple distortions from the cubic geometry: (i) counter- 
rotations of adjacent oxygen octahedra about [111], and 
(ii) relative ionic displacements along [111]. Freezing in 
the counterrotations alone results in a nonpolar insula- 
tor with i?3c symmetry, a structural phase of BiFeOs 
that, by analogy to LiNbOs, is a possible candidate for 
its high-temperature paraelectric phase. Freezing in only 
the polar mode produces a ferroelectric with R3m sym- 
metry. 



B. Electronic and magnetic properties 

In this section, we present the electronic structure of 
BiFeOs computed in the ferroelectric i?3c structure. We 
also briefly discuss results obtained for the cubic per- 
ovskite and centrosymmetric i?3c structures; both are 
possible paraelectric high-temperature phases. 

The single-particle density of states (DOS) for a single 
spin channel and the DOS of the local Fe d states for both 
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FIG. 2: Calculated densities of states (DOS) for R3c BiFe0 3 . 
Upper panels (a) and (b) show the total DOS for one spin 
channel and the local Fe DOS for both spin channels (minor- 
ity spin states are shown as negative) calculated within the 
LSDA, respectively. Lower panels (c) and (d) show the local 
Fe DOS obtained from the LSDA+U method with U eS =2 eV 
and U e ff=4 eV, respectively. The zero is set to the valence 
band maximum. 



spin channels, calculated within the LSDA, are shown in 
Fig. 01 an d b for BiFeC>3 in the R3c structure. Both spin 
channels exhibit identical total DOS, as required for an 
AFM. The structure is insulating, with a small gap of 
0.4 eV in the LSDA calculation. This gap is significantly 
enhanced after application of the LSDA+U method, as 
seen from Fig. [3:; for the small value U e Q=2 eV, the gap 
is 1.3 eV; for C/ C ff=4 eV (Fig. 0i) it further increases to 
1.9 eV. The narrow bands around the Fermi energy arise 
predominantly from the Fe 3c? states (with some oxygen 
hybridization) and are divided into t^g and e g manifolds, 
as expected from their octahedral coordination. Lying 
below and hybridized with the Fe states is the broad pre- 
dominantly oxygen 2p valence band, which also contains 
a significant amount of Bi 6p character. The lowest band 
shown, at ~ 9.5 eV below the valence band maximum, 
is the Bi 6s band. The large amount of occupied Bi 6p 
states is consistent with the picture of the Bi lone pair as 
the driving force of the ferroelectric distortion^ in this 
class of materials. 

Fig. |3 shows the dispersion of the bands in the energy 
range around the gap. As can be seen in the figure, the 
band gap is indirect for LSDA (U c s=0 eV) and U e s=2 eV, 




FIG. 3: Band structure of R3c BiFe03 in the energy region of 
the gap for U cS =0 eV (LSDA, left panel), U eS =2 eV (middle 
panel), and [7 e ff=4 eV (right panel). The valence and conduc- 
tion band edges are indicated by the dashed horizontal lines. 
The high-symmetry k-points are labeled according to Ref.l39l 



with the bottom of the conduction band located at the 
point Z in the rhombohedral Brillouin zone and the top 
of the valence band between T and Z. For a higher value 
of Z7eff=4 eV, the gap remains indirect but the top of 
the valence band shifts to a location between T and F. 
Photoemission spectroscopy would be of use to determine 
the value of U e g best describing the electronic structure 
ofBiFe0 3 . 

The LSDA local magnetic moment at the Fe atoms 
(integrated within a sphere of radius 1.4 A) is ~3.3 fiB, 
in reasonable agr eement with the experimental value of 
3.75 fin (Ref. I4(IL It is reduced from the formal value of 
5 Hb for high spin Fe 3+ because of the finite bandwidth 
of the 3<i states. Use of the LSDA+U method improves 
the agreement with the experiment, enhancing the Fe 
magnetic moment to a value of 3.8 fj,B for U e g=2 eV and 
4.0 hb for E/cff=4 eV. 

For comparison and to discuss the origin of the ferro- 
electric instability, we provide corresponding results for 
the high symmetry cubic perovskite structure. A primi- 
tive lattice constant of 3.87 A is used for all calculations 
and corresponds to a volume per formula unit equal to 
that of the bulk R3c phase; the unit cell is doubled along 
[111] to accommodate the opposite spin polarization of 
the two Fe ions. The single-particle total density of states 
(DOS), calculated within the LSDA, is shown for one spin 
channel in Fig. along with the local density of states 
for both spin channels on one of the Fe ions (Fig. 0Jd). 
Strikingly, BiFeOa is metallic in this structure within the 
LSDA, with the Fermi energy cutting through the nar- 
row band composed of up-spin e g states and down-spin 
tig states. 

As discussed in Sec. |n] the LSDA often fails to de- 
scribe such narrow d bands correctly. Fig.^J: and d show 
the site-projected local DOS for Fe obtained within the 
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FIG. 4: Calculated densities of states (DOS) for cubic BiFe03. 
Upper panels (a) and (b) show the total DOS for one spin 
channel and the local Fe DOS for both spin channels respec- 
tively (minority spin states are shown with a negative sign), 
both calculated within the LSDA. Lower panels (c) and (d) 
show the local Fe DOS obtained using the LSDA+U method, 
with U c fi values of 2 eV and 5 eV. The zero is set to the Fermi 
energy (a,b) or valence band maximum (c,d). 



LSDA+U method for two different values of U c g. Includ- 
ing an effective Coulomb repulsion parameter pushes the 
Fe majority d-bands down and the Fe minority (i-bands 
up in energy, leading to a gap of ~0.5 eV for a relatively 
small U e ff of 2 eV. The fully-gapped insulating phase is 
accompanied by a marked increase in the local Fe mag- 
netic moment from ~1.9 /^b to ^3.9 /ib (indicating a low- 
to high-spin transition). Typical values for U e ff in other 
Fe compounds range between f7 e ff=3 eV and 5 eVj 41 i 42 
the fact that such a small U s opens a gap supports the 
notion that inadequate treatment of strong correlation in 
LSDA results in a metallic state for cubic BiFeOs . 

Finally, we remark that the centrosymmetric i?,3c 
structure, mentioned briefly in Section IIII Al as a pos- 
sible candidate paraelectric phase, is already insulating 
within the LSDA, albeit with a small gap whose mag- 
nitude is found to be extremely sensitive to input pa- 
rameters. As with the cubic perovskite phase, a stable 
insulating solution can be obtained even for small values 
of J7 e fF- We perform a full structural relaxation imposing 
R3c symmetry within the LSDA, and the structural pa- 
rameters resulting from this relaxation are summarized 



in Table IIIII Briefly, compared with i?3c, this arrange- 
ment is found to possess a smaller volume and a larger 
rhombohedral angle. The environment of the A site is 
drastically changed relative to cubic perovskite: each Bi 
site has 3 oxygen neighbors at 2.30 A, 6 at 2.71 A, and 
3 more at 3.21 A. However the local octahedral environ- 
ment of the B site is preserved: each Fe site is octahe- 
drally coordinated with 6 O neighbors at a distance of 
1.97 A. 

TABLE III: Calculated structural parameters (lattice con- 
stant drh; rhomobohedral angle a, and volume fi) of BiFeOs in 
the space group R3c (point group Dsd); the Wyckoff positions, 
referenced to the rhombohedral system, are (2b) Bi(~, ~, ~), 
(2a) Fe(0,0,0), and (6e) 0(x, ± - x,\). 



LSDA 


X 


0.417 


a lh (A) 


5.35 




61.93 


n (A 3 ) 


113.12 



C. Spontaneous polarization 

1. Estimate using Born effective charges 

As discussed above, recent reports of ferroelectric po- 
larizations in high quality BiFeOs thin filmsS*^ exceed 
previous measurements on bulk samples 13 by an order 
of magnitude. The polarizations measured in thin film 
samples are consistent with the observed large atomic 
distortions; 11 ! 12 but apparently inconsistent with an ear- 
lier study of bulk BiFeOa, 13 which reported a modest po- 
larization along [111] of just 6.1 /xC/cm 2 . 

In an effort to shed light on this issue, we examine 
the spontaneous polarization, P, of BiFeOs m its i?3c 
ground state from first principles. Experimentally, P 
corresponds to half the polarization change as the ap- 
plied field is swept through zero from positive to neg- 
ative. For prototypical perovskite ferroelectrics, it has 
been standard to estimate this value by simply summing 
the product of atomic displacements (from a centrosym- 
metric reference structure) and their corresponding Born 
effective charges (BECs). This estimate corresponds to 
computing the Cartesian components of the polarization 
AP a to linear order in the atomic displacements, i.e. 




where Au^ is the displacement of ion j in Cartesian di- 
rection j3, Zj a B is its Born effective charge tensor, and 
is the unit cell volume. The zero subscript generally 
refers to an insulating centrosymmetric reference struc- 
ture (in this case, either cubic perovskite or i?3c). The 



6 



TABLE IV: Born effective charges (BECs) for displacements 
along [111] for BiFeC>3 in the cubic perovskite, R3c, and R3c 
structures. LSDA+U results are obtained using U c s = 2 eV. 
Since cubic perovskite is metallic within LSDA and the R3c 
structure is nearly so, LSDA values are given only for the R3c 
structure. AP, calculated using Eq. Q along a path from the 
cubic structure to R3c, is given for each set of BECs. 





LSDA 


LSDA+U 




Bi Fe O 


Bi Fe 


cubic Zj 
AP 




6.32 4.55 -3.06 
123.1 /iC/cm 2 


R3c Z* 
AP 




4.92 4.25 -3.06 
101.2 /iC/cm 2 


R3c Z* 
AP 


4.28 3.26 -2.50 
84.2 /iC/cm 2 


4.37 3.49 -2.61 
87.3 /xC/cm 2 



spontaneous polarization is obtained from Eq. Q by 
considering a specific structural (or switching) pathway 
parametrized by the change in atomic displacements con- 
necting a centrosymmetric reference structure and P3c. 

Table llVl summarizes our calculated BECs for the three 
different structures studied in this work. The BECs are 
calculated by finite differences: the ions are displaced 
slightly along [111] from their equilibrium positions, and 
the resulting change in polarization is calculated using 
the Berry-phase method described in Sec. [n] The dis- 
placements are chosen to be small enough to ensure the 
validity of the linear treatment in Eq. l[T]l. (Typical dis- 
placements used here are of the order of 0.01 A.) For 
the LSDA+U calculations we use U e s=2 eV, since this 
value gives an insulating solution for all systems. (Higher 
values of U e s would result only in small quantitative 
changes) . 

The highly anomalous values of the BECs computed 
in the cubic structure, and the reduction of these val- 
ues after freezing in the structural distortions, are in 
ageement with former observations in other ferroelectric 
materialsi^^ In particular, the large values of the Bi 
BECs emphasize the important role of this ion as driving 
force of the ferroelectric distortion^ This is in some con- 
trast to the Li ion in the related compound LiNb03, 38 
whose effective charge is found not to deviate appreciably 
from its nominal value. 

In Table IV we also report the polarization difference 
between the cubic (or i?3c) centrosymmetric structures 
and P3c. Using Eq. JTJ and the R3c LSDA BECs, 
84.2 yuC/cm 2 is obtained along [111]^& a value consis- 
tent with the large displacements and BECs. Yet since 
Eq. is only valid for small displacements and because 
the BECs are found to change considerably along the 
path (compare the LSDA+U results for cubic, R3c, and 
R3c in Table IIV|) , this value is only a rough estimate 
of the polarization in BiFeOa. This can be clearly seen 
from the larger values of 123.1 and 101.2 /xC/cm 2 for 



TABLE V: Minimum value of the polarization lattice corre- 
sponding to the R3c structure, computed with different values 
of Ueff- All values are given modulo the polarization quantum. 



t/eff [eV] 


2 4 


P hiC/cm 2 ] 


6.6 2.6 1.2 



AP, also given in Table |IV| obtained using the cubic 
perovskite and R3c BECs, respectively. The accuracy 
could be improved by breaking the path into shorter seg- 
ments and recomputing the BECs at intermediate points 
on the path. However, it is more efficient to compute the 
polarization directly from the Berry-phase theory, as will 
be described in the next section. 



2. Modern theory of polarization 

An accurate quantitative method for computing the 
polarization to all orders in displacement is the so-called 
Berry phase (or "modern" ) theory of polarizationji2i22i2i 
discussed in Sec.|n] In this approach, the electronic con- 
tribution to the polarization is calculated as a geometric 
phase, formally equivalent to the sum of Wannier cen- 
ters of the occupied bands (where each center is assigned 
a charge e)i2£& Because of the Born- von Karman pe- 
riodic boundary conditions employed here, there is an 
ambiguity in the choice of unit cell and the total polar- 
ization may only be determined up to an integer multiple 
of the polarization quantum eR/f2, where e is the charge 
of the electron, R is a lattice vector in the direction of 
polarization, and il is the volume of the unit celli^^^ 
Thus for a given structure, the theory yields a lattice 
of values. The difference in polarization between two 
structures (e.g., a polar structure and its enantiomorph) 
is therefore determined only up to an integer multiple 
of the polarization quantum. To predict which of the 
values would be obtained in a Sawyer- Tower or CV ex- 
periment, it is also necessary to specify a switching path 
along which the system stays insulating in all intermedi- 
ate structural states^ 

In the case of BiFeOs the polarization lattice of the R3c 
ground state is computed to be (6.6 + n ■ 184.2) /iC/cm 2 
along [111] within the LSDA, where n is an integer 
and |eR/f2|=184.2 /iC/cm 2 is the polarization quantum. 
Similar results are obtained with LSDA+U (see TableEJ), 
though a slight shift in polarization with increasing U e g is 
observed. If the point n=0 is selected, near-perfect agree- 
ment is obtained with the measurement of 6.1 /xC/cm 2 
by Teague et ali& over thirty years ago. However, it is 
initially puzzling that none of the allowed values are at 
all close to the estimate of 84.2 /iC/cm 2 made in the pre- 
vious section, using Eq. Q and based on paths involving 
nearby centrosymmetric structures. 

This discrepancy is immediately resolved by comput- 
ing the allowed values of the polarization difference be- 
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tween a centrosymmetric structure and R3c. In this case 
it turns out that the centrosymmetric R3c and cubic per- 
ovskite structures have a polarization not of zero (modulo 
the polarization quantum) but of half of a polarization 
quantum (92.1+71-184.2) i<C/cm 2 along [111] 41 Thus one 
of the allowed values for the polarization difference be- 
tween i?3c and the centrosymmetric reference structures, 
computed with the Berry-phase theory, is 94.7 /iC/cm 2 
{U c ff — 2 eV), agreeing well with the linear estimate of 
the previous section. That the centrosymmetric struc- 
ture has a non-zero value of polarization may at first 
seem counterintuitive. However, it is readily understood 
using the Wannier function reformulation of the Berry- 
phase expression* 2 " Making the conservative assumption 
that the Wannier functions are centered on atoms with 
multiplicity consistent with formal valences of +3 for Bi 
and Fe and —2 for oxygen; taking the origin at a Bi 
atom; and choosing as the additional basis atoms the Fe 
in the center of the ideal cubic perovskite unit cell and 
the three oxygens at the centers of the faces, the Wan- 
nier sum yields a polarization of 92.1 /LtC/cm 2 along [111] 
for the doubled unit cell (using the i?3c volume). This 
value is exactly half of a polarization quantum, the only 
non-zero value allowed by symmetry for a centrosymmet- 
ric structure^ and also exactly the same value that is 
obtained by the exact calculation of the polarization for 
both centrosymmetric structures (cubic perovskite and 
R3c). Alternate choices of origin and basis could result 
in other points on the polarization lattice, all differing by 
a polarization quantum, but would never produce zero for 
these BECs. 

To connect the allowed polarization difference of 
94.7 iiC/cm 2 to a specific path in structural space, we 
compute the polarization from the Berry-phase theory for 
the endpoints and several intermediate structures along 
an idealized "switching path" connecting the positively- 
oriented + (i?3c) with its enantiomorphic counterpart, 
the negatively-oriented — (i?3c) structure with opposite 
polarization, through the centrosymmetric cubic per- 
ovskite structure. Although the actual displacements 
associated with switching will obviously be more com- 
plex, the atomic positions are assumed to travel smoothly 
from positive to negative orientation through the cu- 
bic arrangement. Also, for simplicity these calculations 
assume a rhombohedral angle a=60° for R3c, which 
slightly shifts the polarization quantum to 185.5 /LtC/cm 2 
and the Berry phase polarizations by a few percent. 
Since the cubic perovskite structure is metallic within 
LSDA, we use the results obtained for U c b = 2 eV. 
The results are shown in Fig. For the specific con- 
tinuous insulating path connecting the — (i?3c) lattice 
point corresponding to —2.3 zxC/cm 2 with the +(i?3c) 
point of 187.8 /zC/cm 2 , the polarization is seen to evolve 
smoothly (but nonlinearly) through the centrosymmet- 
ric cubic structure, which has a nonzero polarization. 
Subtracting the two endpoint values gives a polariza- 
tion change of 190.1 /iC/cm 2 and a predicted sponta- 
neous polarization of half that value, or 95.05 /iC/cm 2 . 
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FIG. 5: Change in polarization P along a path from the orig- 
inal R3c structure through the centrosymmetric cubic struc- 
ture to the inverted — (R3c) structure calculated with the 
LSDA+U method and U e s = 2 eV. The possible values of 
P for fixed distortion differ by multiples of the polarization 
quantum, here 185.6 /iC/cm 2 for a=60°. 



(The slight difference from the value of 94.7 /iC/cm 2 re- 
ported above is due to the assumption a=60°.) This 
value is consistent with measurements on (111) oriented 
thin films j& as well as our estimate with Eq. QJ. 



IV. DISCUSSION: POSSIBILITY OF MULTIPLE 
POLARIZATION PATHS 

In the previous section we observed that if the system 
switches along the path indicated in Fig. (or along one 
that can be continously deformed into it) the modern the- 
ory of polarization predicts a measured spontaneous po- 
larization of ~90-100 //C/cm 2 , depending on the choice of 
U e s- In principle however, an infinite set of polarization 
differences between the +(R3c) and — (i?3c) structures is 
possible, provided that a suitable pathway can be found 
to connect any two endpoints. For example, if a switching 
path could be found taking +(i?3c) to — (R3c) through an 
intermediate structure with zero polarization, the mea- 
sured polarization would be 2.3 /iC/cm 2 (for U c e = 2 eV), 
consistent with reports for bulk samples^ Thus although 
the possibility of the small experimental values result- 
ing from poor sample quality cannot be ruled out, it 
is nonetheless interesting that different reported values 
for polarization (both lower, w 6 /iC/cm 2 , and higher, 
> 150 /xC/cm 2 ) may be explained by different switching 
paths. 

Specific paths connecting the smallest +(R3c) polar- 
ization value with the smallest — (i?3c) value have yet 
to be determined. One approach to finding such a path 
would be to identify a centrosymmetric reference struc- 
ture with zero polarization (modulo a quantum), the path 
then being the atomic displacements from i?3c to this 
structure. To get a zero (or integer-quantum) polariza- 
tion from only a small distortion of the perovskite struc- 
ture would require a somewhat drastic rearrangement of 
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the centers of the Wannier functions. In the present case 
this may be facilitated by the known multiple valences 
possible for Bi and Fe. For example if the Bi ions were 
to acquire an average formal charge of +4 (likely if a dis- 
proportionation to Bi 3+ and Bi 5+ occured, leaving the 
Fe with a +2 charge), the polarization of the cubic phase 
would be zero. Further, by combining a forward switch- 
ing path between i?3c enantiomorphs through an integer- 
polarization structure with another back through a half- 
quantum polarization structure (e.g., i?3c), the BiFeOa 
crystal could be taken to itself with a net transport of 
electrons across the system. This suggests the intriguing 
possibility of an insulating crystal with nonzero electronic 
conductivity. 

In constructing paths corresponding to different po- 
larizations, we do not wish to claim knowledge of the 
actual experimental switching mechanism, which is cer- 
tainly much more complicated. In reality, switching is 
thought to occur via domain wall motion, where the key 
physics is associated with polarization reversal under an 
applied electric field at interfaces between positively- and 
negatively-oriented domain walls4& 

Finally, an alternative explanation for the different 
polarizations reported in various film and bulk samples 
would be that different crystal structures are epitaxially 
stabilized in the films which in turn possess significantly 
different polarizations. An alternative tetragonal struc- 
ture that satisfied this requirement was discussed in a 
previous paper£ To clarify these issues further experi- 
mental work will be necessary to determine the sensitiv- 



ity of the structure of thin film BiFeC>3 to substrate and 
growth orientation. 



V. CONCLUSIONS 

In summary, BiFeC>3 is a material of unusual interest 
both as a potentially useful multiferroic, and with respect 
to its fundamental polarization behavior. A wide range 
of measured polarization values have been reported, all 
of which are apparently permitted within the modern 
theory of polarizationi2i22i2i by the lattice character of 
the polarization. Since some of the observed values of 
polarization can only be explained by switching through 
structures in which the ions change their valence states, 
such behavior, if experimentally verified, might be unique 
to multiferroics, in which the magnetic transition metals 
are able to adopt multiple values for d orbital occupancy. 
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